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ABSTRACT 

We describe how to generate a longitudinal brightness map for a tidally locked exoplanet from its 
phase function light curve. We operate under a number of simplifying assumptions, neglecting limb 
darkening/brightening, star spots, detector ramps, as well as time- variability over a single planetary 
rotation. We develop the transformation from a planetary brightness map to a phase function light 
curve and simplify the expression for the case of an edge-on system. We introduce two models — 
composed of longitudinal slices of uniform brightness, and sinusoidally varying maps, respectively — 
which greatly simplify the transformation from map to light curve. We discuss numerical approaches 
to extracting a longitudinal map from a phase function light curve, explaining how to estimate the 
uncertainty in a computed map and how to choose an appropriate number of fit parameters. We 
demonstrate these techniques on a simulated map and discuss the uses and limitations of longitudinal 
maps. The sinusoidal model provides a better fit to the planet's underlying brightness map, although 
the slice model is more appropriate for light curves which only span a fraction of the planet's orbit. 
Regardless of which model is used, we find that there is a maximum of ~ 5 free parameters which can 
be meaningfully fit based on a full phase function light curve, due to the insensitivity of the latter 
to certain modes of the map. This is sufficient to determine the longitudes of primary equatorial 
hot-spots and cold-spots, as well as the presence of secondary maxima/minima. 
Subject headings: methods: data analysis — (stars:) planetary systems — 



1. INTRODUCTION 

Observations of secondary eclipse s in exoplanetary sys- 
tems, startin g with HP 209458b (iDeming et all I2005T) 
and TrES-1 (ICharbonneau et al ] I2005I K have made it 
possible to measure the integrated day-side flux of hot 
Jupiters (for a review of t ransiting exoplanet science see 
ICharbonneau et "aTJ l2007t l . By carefully studying of the 
shape of the ingress and egress of secondary eclipses, it 
should e ventually be possible to map the day-side of such 
planets (|Williams et al.ll2006l : lRauscher et al.ll2007aD . To 
characterize the planet's longitudinal temperature pro- 
file at all longitudes, however, observations must be 
made at a variety of points in the planet's orbit. The 
firs t successful observation s of this sort we r e rep orted 
bv lHarrington et all (|2006D . iKnutson et ail (|2007t ) and 
iCowan et al.l (120071). The un prec edented quality of the 
data in IKnutson et all (|2007f) and IKnutson et alj 1)20081 ) 
made it possible to model the planet not merely as day 
and night hemispheres, but to divide the planet into lon- 
gitudinal slices, hence producing the first (albeit coarse) 
maps of an exoplanet under restricted assumptions. 

In this paper we elaborate on the inversion techniques 
which one can use to obtain a longitudinal brightness 
map from the light curve of phase variations. Such maps 
promise to be powerful diagnostic tools for simulations 
of hot Jupiter atmospheric dynamics because they are 
nearly model- independent. 

1.1. Simplifying Assumptions 

The detailed study of coupled radiative transfer and 
dynamics in the atmospheres of hot Jupiters is a tremen- 
dously complex science (for a current review of the field, 
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see lShowman et afll2007l) . Certain models of hot Jupiter 
atmospheres predict significant variability in the inte- 
grated brightness of the planets (Rausch er et al.lr2007bl 
and references therein). Although it may be possible to 
glean useful information from the phase function light 
curves of planets with such variable atmospheres, we 
elect in this Letter to ig nore time variability, i n line with 
the detailed modeling of lCooper fc Showm an ( 2005( 1 . We 
also choose to ignore limb darkening/brightening since 
most atmospheric models of hot Jupiters predict infrared 
photon absorption lengths much shorter than the scale 
height of temperature variations. In any case, the pres- 
ence of limb darkening would not significantly change the 
present analysis since the limb does not contribute much 
to the integrated light of the planet (we verify this in 
§ 3.4). 

Star spots can plague phase function observations since 
they change the overall brightness of a planet/star sys- 
tem as they rotate into and out of view. Fortunately, 
star spots vary on longer timescales (~weeks) than the 
orbital periods of currently known transiting exoplanets 
(~days), and have larger variations in the optical/near- 
IR which can be used t o characterize and su btract the 
star spot variation (eg: IKnutson et al.l [2008) . Further- 
more, some instrum ents on the Spitzer Space Telescope 
(|Werner et al.l2004f ) exhibit detector ramps which distort 
phase function observati ons, especially near t he begin- 
ning of a time series (eg: IKnutson et al.l [2007h ■ Both of 
these observational challenges can be overcome through 
longer observations: full orbits or more. For the purposes 
of this Letter, we assume that any — and all — variations 
in the infrared brightness of an extrasolar planetary sys- 
tem are due to changes in the flux from the planet as 
different portions of the planet rotate in and out of view. 
SfMlar variability, which contributes to the photometric 
scatter in the light curve, will likely be the limiting factor 
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in future phase function analysis. 

The mapping technique described here is predicated 
on the known rotation rate of the planet. As such, it is 
only applicable to tidally locked planets, where the rota- 
tion and orbital periods are identical. Most hot Jupiters 
are tidally locked so this is not a problematic restric- 
tion. Eccentric hot Jupiters are thought to be in pseudo- 
synchronous rotation, in which case their rotational pe- 
riod can be derived from their orbital parameters. It is 
unlikely that eccentric planets have steady-state atmo- 
spheric dynamics and hence, although the light curve in- 
version methodology described in this Letter may be ap- 
plicable, it is not clear that the result would be a "map" . 

2. FROM MAPS TO LIGHT CURVES 

The object of this Letter is to demonstrate how to 
transform an observed light curve, F(£), into a longi- 
tudinal brightness map of the planet, where the planet's 
phase, £, is the angle — in the plane of the planet's orbit — 
between the planet, its host star and the planet's position 
at superior conjunction (£ = at superior conjunction; 
£ = 7r at inferior conjunction). For edge-on systems, £ 
corresponds to the observer-planet-star angle (£ = at 
secondary eclipse; £ = tt at transit). In the interest of 
simplicity, we ignore the dips in light due to transit and 
secondary eclipse which occur in transiting systems, al- 
though these are crucial in pinning down the absolute 
flux of the planet (for non-transiting planets F(£) is only 
known to within an additive constant). 

If a and (3 are the longitude and latitude on the plan- 
etary disk as seen by an observer and I is the intensity 
map of the planet, the total flux from the planet as seen 
by that observer is: 

F(£)= / I(a, (3) cos a sin 2 (3dadf3. (1) 

For a tidally locked planet, it is possible to define longi- 
tude, <p, and latitude, 9, in the planet's rotating frame, 
such that <j) = at the sub-stellar point, 9 — at the 
planet's north pole, and <p increases in the direction of 
rotation of the planet. The rotation that relates (a, (3) to 
(4>, 9) can be expressed in terms of Euler angles: <f> = tt/2, 
6 = 7r/2 — i, "J" = £ — 7r/2, where i is the inclination of the 
planet's orbit (i = for a face-on orbit; i = f for an edge- 
on orbit). The rotation leads to a system of three coupled 
equations which can be solved for a(4>, 9) and /?(</>, 9). By 
inserting these expressions into Equation [TJ one obtains 
an equation for the observed flux from the planet as a 
function of the specific intensity of the planet at differ- 
ent longitudes and latitudes — in effect, a transformation 
from a brightness map to a light curve: 1(9, 0) — ► F(£). 

2.1. Edge-on Orbits 

Although the formalism above is sufficient to produce 
a light curve from a two-dimensional planetary map, it 
is instructive to further constrain the problem by consid- 
ering planets with inclinations of precisely tt/2, with the 
understanding that the resulting solutions should be ap- 
proximately correct for transiting planets. In the worst- 
case scenario of a transiting system with i = 80° and 
7(0, 9) constant in 9, the edge-on approximation leads to 
5% errors in F(£). For edge-on orbits, the rotation re- 
lating (a, (3) to ((f), 9) simplifies to a = + £, and (3 = 9. 



The expression for F(£) can then be written as 

/"TT r- f-Hr/2 

F(£) = / / 7(0, 9) cos(0 + sin 2 9d<j)d9. (2) 

JO J-t-Tr/2 

There are no current observations which can constrain 
the ^-dependence of /, but for edge-on orbits the lati- 
tudinal dependence of the intensity is unimportant since 
one can define J{4>) = Jq I(4>, 9) sin 2 9d9, which repre- 
sents the flux contribution from an infinitesimal slice of 
the planet when viewed face-on. The integrated flux from 
the planet at a given point in its orbit is then 

F(0= / C+V2 J(0)cos(0 + £)d0. (3) 

It is useful to think of this integral as a convolution 

= J J{(j))K((j), £)d0, with the piece-wise defined 
kernel K: 

K(rh cos ^ + if cos ^ + - U\ 

A ^'^"\0 if cos(0 + f) < 0. W 

The kernel, which represents the response of the phase 
function to a delta function in J, is very broad, with a 
full width at half-maximum of 27r/3. 

3. FROM LIGHT CURVES TO MAPS 

In the previous section we developed an analytic ex- 
pression, Equation [3l for the convolution J(0) — > F(£). 
There is no closed expression for the deconvolution, 
— * J((f>), and there is no guarantee that solutions 
are unique. Furthermore, two problems arise in prac- 
tice: the light curve is only sampled at discrete values 
of £ which may not span a full planetary rotation, and 
the measurements of i r (£) are not arbitrarily precise but 
instead have associated uncertainties o~p. Given these 
realities, it is useful to develop model maps which sim- 
plify the integral in Equation [3] and then use numerical 
methods to solve the problem or — better yet — allow for 
direct inversion. 

3.1. N-Slice (AKA "Orange Slice") Model 

The planet is divided into N longitudinal slices of 
width Ac/) = 2ir/N. Each slice has a uniform intensity in 
both longitude and latitude. This flux distribution is not 
ruled out by the observations and — more importantly — 
smoothing the steps does not significantly change the 
light curve, provided the total flux from each slice and 
their brightness-weighted longitude are unchanged. The 
jth slice is centered on <f>j = 4>q + jAcf), where the phase 
offset, < 4>q < A</>, is useful to accommodate slight 
discrepancies between the light curve maximum and su- 
perior conjunction. The intensity map for the planet is 
given by: 

{J if 0o - Acj)/2 <cj)<(j)o + A<j)/2 
Ji if 0i - A(j)/2 < (j) < 0i + A0/2 
: 
Jat_i if w _i - A0/2 < < 0at_i + A0/2. 

(5) 

Since, in practice, one is only ever concerned with com- 
paring the model phase function to data at a finite num- 
ber of discrete phases, the transformation from N-slice 
map to light curve can be expressed in matrix form: 
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F = GJ . The matrix G is defined as G = sin a+ — sin a_ , 
where a+ and a_ represent the leading and trailing edges 
of the slices: 



a± 



cos 



(max[cos(& + ± A</>/2), 0]) . (6) 



3.2. Sinusoidal Model 

Sinusoidal basis maps have the advantage of producing 
sinusoidal light curves of the same frequency and phase 
offset. If a planet map is composed of sinusoids, J(4>) = 
A-o +YljLi Aj cos(j(f)) + Bj sin(j(/)), the light curve is sim- 
ply given by F(£) = F a + £f =1 Cj cos(j'0 + D s sin(j£). 
The coefficients of F are related to those of J by: 
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where j must be even. Sinusoidal modes with odd j 
(other than j = 1) do not have a phase function signa- 
ture. In Figure [1] we show the light curve contributions 
for a handful of sinusoidal modes, assuming that all of 
the modes have the same amplitude in J. The higher fre- 
quency modes are strongly suppressed due to the broad 
smoothing kernel. This low-pass filter limits the number 
of modes which can be meaningfully fit with a given light 
curve. The uncertainty in a sinusoidal mode in the light 
curve is related to the uncertainty in the map sinusoidal 
modes by Equation [7] (eg: ac t — {n/2)aA 1 ). 




Fig. 1. — The phase function for different sinusoidal maps. Odd 
sinusoidal maps (starting with J(<j>) = cos(3</>)) are invisible due to 
symmetry. 



3.3. Inversion Techniques 

Since both the N-slice and sinusoidal models described 
above provide computationally efficient ways to gener- 
ate light curves from maps, one can use a \ 2 fitting rou- 
tine (Markov Chain Monte Carlo, Levenberg-Marquardt, 
etc.) to produce a map from a given light curve. It is sim- 
ply necessary to demand that J(<fi) be strictly positive. 
These techniques have the advantage of naturally pro- 
ducing error estimates for the resulting map. Although 



a unique best-fit map can always be determined in this 
way, the uncertainty in the fit parameters may be very 
large if the number of parameters is not commensurate 
with the signal-to-noise of the light curve. 

For the N-slice model, the uncertainty balloons for all 
the Ji if too many free parameters are used. We there- 
fore suggest running multiple fits with different numbers 
of slices. When the addition of a slice does not improve 
the x 2 , one has achieved the best model that the data can 
support. Instead of repeating the fit with fewer slices, one 
can apply smoothing to the map in the form of a Bayesian 
prior (eg: iKnutson et alJ 120071 ) . For well-designed pri- 
ors, however, there is no fundamental difference between 
models with large N and long smoothing lengths, ver- 
sus models with smaller TV and little or no smoothing. 
In the interest of simplicity we recommend using fewer 
slices rather than smoothing. 

An observed light curve can be quickly deconvolved 
into sinusoidal maps by determining its Fourier compo- 
nents, Fq, C and D, then converting them via the Equa- 
tion [7J It is expedient to assume that there is no power 
in the odd modes (other than j = 1) to avoid degenerate 
solutions, but this may lead to a systematic error in the 
model map, depending on how much power is present 
in these modes in the real map. The uncertainty in the 
map parameters may be determined using a x 2 routine 
or Monte Carlo analysis. Since the sinusoidal model has 
linearly independent modes, only the uncertainty in the 
highest-frequency modes explodes when too many modes 
are considered. As a rule of thumb, one should truncate 
the Fourier series once the uncertainty in coefficients be- 
comes greater than the coefficients themselves. 

3.4. Inversion Example 

We now turn to an example map and test the abil- 
ity of the algorithms described above to recover the cor- 
rect features of this map. The top panel of Figure [2] 
shows the brightness map computed from a snapshot o f 
the atmospheric dynamics model of ICho et aTT ([2003) . 
We performed an analogous test — and obtained com- 
parable results— using a snapshot from the model of 
I Showman et al.l ((2008) , which we do not include here 
in the interest of space. The brightness map was gen- 
erated by treating each pixel of the temperature map 
as a blackbody and computing an associated intensity 
at 8^m. The bottom panel of Figure [5] shows the inte- 
grated longitudinal brightness map, J ((/)). Note that the 
sin 9 term in the integral for J attenuates the flux con- 
tribution from the poles of the planet. Also shown in the 
bottom panel are the best fit maps for the N-slice and 
sinusoidal models. The light curve associated with the 
map, as well as the best-fit light curves for the models, 
are shown in Figure [3] The map J was converted to an 
idealized light curve, and mock observations (comprised 
of 100 data points) were generated by removing the seg- 
ments of the light curve corresponding to the transit and 
secondary eclipse, then scaling the planet/star flux ra- 
tio and the photometric uncer tainties to roughly match 
those of IKnutson et alJ (|2007f ). Both models reproduce 
the features of the map, as well as the light-curve. The 
best-fit 5-slice model was determined using a Levenberg- 
Marquardt ^-minimization routine. The sinusoidal map 
was determined by decomposing the light curve into si- 
nusoidal components to j = 2, then converting the coef- 
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ficients using Equation [71 The uncertainties — estimated 
by Monte Carlo analysis — in the j = 4 terms are larger 
than their amplitudes and are therefore ignored. The 
insensitivity of phase functions to the first and most im- 
portant odd mode, j = 3, leads to ~ 10% errors in the 
resu lting map, based on t he maps of ICho et al.1 ()2003l ) 
and I Showman et alj (|2008[ ). We model the effect of limb 
darkening by adding [0.8 + 0.2 sin 6] to the expression for 
J{(f)) and adding [0.8 + 0.2cos(<^) + £)] to the integrand 
of Equation [3] We find the resulting light curve to differ 
by less than 0.2%, justifying our decision to neglect limb 
darkening in the formalism above. 




Fig. 2. — The top panel shows a brightness map (darker is 
cooler) based on the atmospheric dynamics model of ICho et all 
(t2003T) . with T min = 1020 K and T max = 1269 K. In the lower 
panel, the solid line is J(<f>), the histogram and associated error 
bars represents the best-fit 5-slice model, and the gray band is the 
lcr confidence interval for the sinusoidal map. Note that positive 
<f> are to the left of the plot, to facilitate comparison with the light 
curve in Figure [3] 




n/2 



Tt/2 



Fig. 3. — The phase function light curve for the planet map 
shown in Figure [2] convolved with photometric scatter comparable 
d. ( 2 



4. DISCUSSION & CONCLUSIONS 

The best current light curves can be represented by the 
function = F Q + C\ cos f + £>i sin £ + C 2 cos(2£) + 
£>2 sin(2£) which can be directly translated into a lon- 
gitudinal brightness map using Equation [7J The j = 3 
modes should cancel out by symmetry so their presence 
in the light curve would indicate systematic errors; j > 4 
modes are generally lost in the noise. The longitudes 
and amplitudes of the primary hot-spots and cold-spots 
can be determined from the j = 1 terms, while the rel- 
ative strength of the j = 1 and j = 2 modes indicates 
whether there are secondary local maxima/minima. By 
the same token, a 4-slice model with variable phase offset 
should be sufficient to model most phase function light 
curves. For light curves with incomplete phase coverage, 
the uncertainty in the sinusoidal map is the same at all 
longitudes whereas the N-slice model naturally has larger 
uncertainties for slices which were visible for less time, 
an intuitive and desirable property. So far only half- 
orbits have been allocated to phase function studies, but 
the warm Spitzer mission will provide a perfect opportu- 
nity to obtain light c urves spanning full planetary orbits 
(|Deming et al.ll2007f) . 

The contribution of sinusoidal modes to the observed 
light curve decreases precipitously with j. Not only does 
the J(<p) — > F^) transformation suppress the terms as 
l/(j 2 — 1), but the intrinsic power of these modes in the 
underlying m ap might drops as ~ j~ 2 (as is the case with 
themap from ICho et al.l [20031 ). As a result, the j = 4 
modes in light curves will be ~ 20 times weaker than the 
j = 2 modes and will likely remain undetectable with 
JWST. Even when higher quality light curves are even- 
tually obtained — and assuming that stellar variability is 
not the limiting factor — the physical significance of the 
model map would be questionable due to the insensitiv- 
ity of the phase function to odd sinusoidal modes. The 
power in these modes can be constrained theoretically 
by dynamical atmospheric models and observationally 
through secondary eclipse mapping, which promises to 
be feasible with JWST. 



to lKnut son et al. ( 2007). The solid and dashed lines shows the light 
curves of the best-fit N-slice and sinusoidal models, respectively. 
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